rm(list = ls())
load("dataBJPOLS.RData")
last.cases[, dups := .N, by = c("this_id")]

last.cases$type <- NA
last.cases$type[last.cases$dups == 1 & last.cases$time_period == 1] <- 1
last.cases$type[last.cases$dups == 1 & last.cases$time_period == 2] <- 2
last.cases$type[last.cases$dups == 2] <- 3

last.cases <- last.cases[last.cases$time_period == 1, ]

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + as.factor(noteli) + regis_before + vote_pre"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop)"

## Mediation
outc.1 <- "I(vote_2016 * regis_before)"
med.1 <- "vote_post"

form.1m <- formula(paste(med.1, "~", endo.1, "+", time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
form.2m <- formula(paste(outc.1, "~", med.1, "+", endo.1, "+", time.controls, "+", demo.controls, "+", case.controls, "|", endo.1, "+", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

results <- list()
boot.run <- T

if(boot.run) {
  set.seed(1231)
  for(i in 1:500) {
    out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
    last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
    
    m1a3m1 <- ivreg(form.1m, data = last.cases.boot)
    m1a3m2 <- ivreg(form.2m, data = last.cases.boot)
    m1a3m3 <- ivreg(form.3, data = last.cases.boot)
    
    direct.effect <- coef(m1a3m2)["pti"]
    total.effect <- coef(m1a3m3)["pti"]
    indirect.effect <- total.effect - direct.effect
    
    results[[i]] <- c(direct.effect, indirect.effect, total.effect)
  }
  
  SEs <- round(apply(do.call('rbind', results), 2, sd), 2)
}

m1a3m1 <- ivreg(form.1m, data = last.cases)
m1a3m2 <- ivreg(form.2m, data = last.cases)
m1a3m3 <- ivreg(form.3, data = last.cases)

coef(m1a3m2)["pti"] + coef(m1a3m1)["pti"] * coef(m1a3m2)["vote_post"]
direct.effect <- coef(m1a3m2)["pti"]
total.effect <- coef(m1a3m3)["pti"]
indirect.effect <- total.effect - direct.effect

main <- c(direct.effect, SEs[1], indirect.effect, SEs[2], total.effect, SEs[3])
round(main, 3)



